scaled_coef_cluster <- function(ModelResults, data, clusterid, subvar = FALSE, vars = NULL) {
       library(multiwayvcov)
       library(lmtest)
       library(dplyr)
       require(ggplot2)
       library(gridExtra)
       library(tidyr)
       library(stringr)
       library(purrr)
       
       cluster <- data[, clusterid]
       
       vcov_cluster = lapply(ModelResults, function(x) FUN = cluster.vcov(x, 
                                                                          cluster))
       modcoef = Map(function(x, y) coeftest(x, y), x = ModelResults, y = vcov_cluster)
       modcoef <- lapply(modcoef, function(x) FUN =round(x, digits = 4))  
       modelcoef = lapply(modcoef, function(x) FUN = x[, "Estimate"])
       modelse = lapply(modcoef, function(x) FUN = x[, "Std. Error"])
       varnames = lapply(ModelResults, function(x) FUN = names(x$coefficients))
       dfplot <- list()
       for (i in 1:length(modelcoef)){
              ad = 0.5 / modelse[[i]]
              modelcoef[[i]] = modelcoef[[i]]* ad
              modelse[[i]] = modelse[[i]]*ad 
              
              ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
              yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
              ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
              yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
              dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                                       ylo90, yhi90,
                                       Model = paste('Model',seq_along(modelcoef)[i]))
       }
       
       coefficientdata = do.call(rbind,dfplot)
       colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                                      "Low95CI", "High95CI", 
                                      "Low90CI", "High90CI",
                                      "Model.Name")
       df <- coefficientdata %>%
              dplyr::mutate(Variables = as.character(Variables)) %>%
              dplyr::arrange(Variables)
       
       if(subvar == TRUE){
              df <- df[grep(paste(vars, collapse = "|"), df$Variables),]
              #re-order the variables: if factor, the order of factor levels is used; if character, an alphabetical order ist used
              df$Variables <- factor(df$Variables, levels = vars)
       } 
       
       p = ggplot(df, aes(colour = Model.Name)) + 
              geom_hline(yintercept = 0, lty = 2) +
              geom_linerange(aes(x = Variables, ymin = Low90CI,
                                 ymax = High90CI),
                             lwd = 2, position = position_dodge(width = 1/2)) +
              geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                                  ymax = High95CI, shape = Model.Name),
                              lwd = 1, position = position_dodge(width = 1/2),
                              fill = "WHITE") +
              coord_flip() + 
              theme_bw() + xlab("") + ylab("Rescaled Coefficient") + 
              theme(legend.position="bottom",
                    legend.title=element_blank(),
                    axis.text = element_text(size=14),
                    text = element_text(size=14),
                    plot.title = element_text(hjust = .5, size = 16, face = "bold"))
       
       return(p)   
} 


scaled_coef_cluster2 <- function(ModelResults, data, clusterid, subvar = FALSE, vars = NULL) {
        library(multiwayvcov)
        library(lmtest)
        library(dplyr)
        require(ggplot2)
        library(gridExtra)
        library(tidyr)
        library(stringr)
        library(purrr)
        
        cluster <- data[, clusterid]
        
        vcov_cluster = lapply(ModelResults, function(x) FUN = cluster.vcov(x, 
                                                                           cluster))
        modcoef = Map(function(x, y) coeftest(x, y), x = ModelResults, y = vcov_cluster)
        modcoef <- lapply(modcoef, function(x) FUN =round(x, digits = 4))  
        modelcoef = lapply(modcoef, function(x) FUN = x[, "Estimate"])
        modelse = lapply(modcoef, function(x) FUN = x[, "Std. Error"])
        varnames = lapply(ModelResults, function(x) FUN = names(x$coefficients))
        dfplot <- list()
        for (i in 1:length(modelcoef)){
                ad = 0.5 / modelse[[i]]
                modelcoef[[i]] = modelcoef[[i]]* ad
                modelse[[i]] = modelse[[i]]*ad 
                
                ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
                yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
                ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
                yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
                dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                                         ylo90, yhi90,
                                         Model = paste('Model',seq_along(modelcoef)[i]))
        }
        
        coefficientdata = do.call(rbind,dfplot)
        colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                                       "Low95CI", "High95CI", 
                                       "Low90CI", "High90CI",
                                       "Model.Name")
        df <- coefficientdata %>%
                dplyr::mutate(Variables = as.character(Variables)) %>%
                dplyr::arrange(Variables)
        
        if(subvar == TRUE){
                df <- df[grep(paste(vars, collapse = "|"), df$Variables),]
                #re-order the variables: if factor, the order of factor levels is used; if character, an alphabetical order ist used
                df$Variables <- factor(df$Variables, levels = vars)
        } 
        
        p = ggplot(df, aes(colour = Model.Name)) + 
                geom_hline(yintercept = 0, lty = 2, color = "#FFFFFF") +
                geom_linerange(aes(x = Variables, ymin = Low90CI,
                                   ymax = High90CI),
                               lwd = 2, position = position_dodge(width = 1/2)) +
                geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                                    ymax = High95CI, shape = Model.Name),
                                lwd = 1, position = position_dodge(width = 1/2),
                                fill = "WHITE") +
                coord_flip() + 
                theme_bw() + xlab("") + ylab("Rescaled Coefficient") + 
                theme(legend.position="bottom",
                      legend.title=element_blank(),
                      axis.text = element_text(size=14),
                      text = element_text(size=14),
                      plot.title = element_text(hjust = .5, size = 16, face = "bold"))
        
        return(p)   
} 


theme_black2 = function(base_size = 18, base_family = "") {
        
        theme_grey(base_size = base_size, base_family = base_family) %+replace%
                
                theme(line = element_blank(),
                      # Specify axis options
                      axis.line = element_blank(),  
                      axis.text.x = element_text(size = base_size*0.8, color = "#FFFFFF", lineheight = 0.9),  
                      axis.text.y = element_text(size = base_size*0.8, color = "#FFFFFF", lineheight = 0.9),  
                      # axis.ticks = element_line(color = "#FFFFFF", size  =  0.2),  
                      #axis.text = element_blank(),
                      # axis.ticks = element_blank(),
                      title = element_text(color = "#FFFFFF"),
                      axis.title.x = element_text(size = base_size*0.8, color = "#FFFFFF", margin = margin(0, 10, 0, 0)),  #t = 0, r = 0, b = 0, l = 0
                      axis.title.y = element_text(size = base_size*0.8, color = "#FFFFFF", angle = 90, margin = margin(0, 10, 0, 0)),  
                      axis.ticks.length = unit(0.3, "lines"),   
                      # Specify legend options
                      legend.title = element_blank(),
                      legend.background = element_rect(color = NA, fill = "#2B2E34"),  
                      legend.key = element_rect(color = "#FFFFFF",  fill = "#2B2E34"),  
                      legend.key.size = unit(1.2, "lines"),  
                      legend.key.height = NULL,  
                      legend.key.width = NULL,      
                      legend.text = element_text(size = base_size*0.8, color = "#FFFFFF"),  
                     # legend.title = element_text(size = base_size*0.8, face = "bold", hjust = 0, color = "#FFFFFF"),  
                      legend.position = "bottom",  
                      legend.text.align = NULL,  
                      legend.title.align = NULL,  
                     # legend.direction = "vertical",  
                      legend.box = NULL, 
                      rect = element_blank(),
                      # Specify panel options
                      panel.background = element_rect(fill = "#2B2E34", color  =  NA),  
                      # panel.border = element_rect(fill = NA, color = "#FFFFFF"),  
                      #panel.grid.major = element_line(color = "grey35"),  
                      # panel.grid.minor = element_line(color = "grey20"),  
                      panel.spacing = unit(0, "lines"),   
                      # Specify facetting options
                      # strip.background = element_rect(fill = "grey30", color = "grey10"),  
                      strip.text.x = element_text(size = base_size*0.8, color = "#FFFFFF"),  
                      strip.text.y = element_text(size = base_size*0.8, color = "#FFFFFF",angle = -90),  
                      # Specify plot options
                      plot.background = element_rect(color = "#2B2E34", fill = "#2B2E34"),  
                      plot.title = element_text(size = base_size*1.2, color = "#FFFFFF"),  
                      plot.margin = unit(c(0, 1, 1, 0), "lines")
                      
                )
        
}




